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ABSTRACT 

Infrared radiation emitted from an asteroid surface causes a torque that can signifi¬ 
cantly affect rotational state of the asteroid. The influence of small topographic fea¬ 
tures on this phenomenon, called the YORP effect, seems to be of utmost importance. 
In this work, we show that a lateral heat diffusion in boulders of suitable sizes leads 
to an emergence of a local YORP effect which magnitude is comparable to the YORP 
effect due to the global shape. We solve a three-dimensional heat diffusion equation in 
a boulder and its surroundings by the finite element method, using the FreeFem-|--|- 
code. The contribution to the total torque is inferred from the computed tempera¬ 
ture distribution. Our general approach allows us to compute the torque induced by 
a realistic irregular boulder. For an idealized boulder, our result is consistent with an 
existing one-dimensional model. We also estimated (and extrapolated) a size distribu¬ 
tion of boulders on (25143) Itokawa from close-up images of its surface. We realized 
that topographic features on Itokawa can potentially induce a torque corresponding 
to a rotational acceleration of the order 10“^radday“^ and can therefore explain the 
observed phase shift in light curves. 

Key words: minor planets, asteroid: individual: (25143) Itokawa - methods: numer¬ 
ical. 


1 INTRODUCTION 

The Yarkovsky-O’Keefe-Radzievskii-Paddack effect is the 
torque caused by the infrared emission from an asteroidal 
surface that has an influence on a rotational state of an as¬ 
teroid ( Rubincam||2000 1. It is now widely recognized as an 
important factor, affecting the evolution of rotational states 
of asteroids alongside mutual collisions and tidal torques. 
The YORP effect helped to explain numerous observed phe¬ 
nomena, such as a spin axis alignment of asteroids in the Ko- 


ronis family | Vokrouhlicky et al. 20031, a non-maxwellian ro¬ 


tational frequency distributions of small main-belt asteroids 
| |Pravec et al.|200^ or significant binary asteroid population 
among near-Earth objects ( Walsh et al.||20i2 K 

Even a direct evidence of a non-gravitational torque has 
been found. A phase shift in light curves has been measured 
for a few asteroids that cannot be explained by a solely grav¬ 
itational model — (1862) Apollo ||Kaasalainen et al.||2007l, 
(54509) YO RP dLowry et al.|2007||Taylor et al.|2007|), (1620) 
Geographos ( Durech eT^irpOOS I, (3103) Eger Durech et al. 
2012), and finally, (25143) Itokawa (Lowry et al. 2014). 


There is no known mechanism except the YORP effect that 
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could explain the quadratic trend in rotational phase ob¬ 
served for these asteroids. 

The asteroid Itokawa has been a suitable candidate 
for a detection of YORP effect for its highly asymmetric 
shape and its favourable position among near-Earth objects. 


Vokrouhlicky et al. (2004 1 predicted a measurable accelera¬ 
tion of rotation of the order of 10“^ rad day based on the 
shape model derived by radar ranging. Itokawa was then 
a target of the Hayabusa spacecraft in 2005 and a state- 
of-the-art shape model of the asteroid was constructed from 


silhouette images (Gaskell et al. 20061. The torque computed 


using the latter model would lead to a significant deceler - 
ation —(1.8 to 3.3) x 10“^radday“^ ( Sch^res et al.|[2007 |. 
However, the measured phase shift in light cur ves revealed 
an acceleration -|-(0.35±0.04) x 10“^ radday“^ |Lowry et al. 


2014). Theoretical models did not predict even the sign of 


the effect correctly. This discrepancy between the observed 
and predicted change of the angular frequency is concerning 
and only one viable explanation has been put forward to 
date. 

The observed rotational acceleration could be at¬ 
tributed to density inhomogeneities in the asteroid. [Scheere^ 
and Gaskell 120081 showed that the YORP effect on Itokawa 


is indeed sensitive to the position of the center of mass. 
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2 P. Sevecek et al. 


Based on the measured acceleration, Lowry et al. (20141 
computed the required offset between the center of mass 
and the center of figure to be ~21m. Such offset indi¬ 
cates that the asteroid might consist of two parts with 
substantially different densities — (2850 ± 500) kg m“^ and 

(1750 ± 110) kgm-^ _ 

The deceleration predicted by Scheeres et al. (20071 was 
computed from the shape model with ~ 50000 facets. Cal¬ 
culations of the effect with a more detailed shape lead to an 
even bigger deceleration. According to Breiter et al. (20091, 
the deceleration does not show any sign of convergence with 
increasing resolution, implying that even sub-meter sized 
surface features possibly have a non-negligible influence. For 
the highest available shape resolution, their model predicted 
the deceleration —5.5 x 10“^radday“^. 


Lowry et al. (20141 employed the advanced thermophys¬ 
ical model ( Rozitis and Green||2011 2012 20131, including 
the effect of thermal infrared beaming and the global self¬ 
heating of the asteroid. By varying the distribution of rough 
surface in patchy ways, the model showed a change of angu¬ 
lar velocity (—1.80±1.96) x 10“^radday“^ (see Fig. 5 in the 
cited paper). Remarkably, an acceleration can be achieved 
even without a shift of the center of mass. However, this 
result was obtained only in 16.5% cases and the roughness 
distribution corresponding to these cases seems rather arti¬ 
ficial. 

There is a problem that shapes with surface features of 
sub-meter sizes cannot be easily included in existing mod¬ 
els of the YORP effect. There are several reasons for this 
limitation. First, numerical YORP models typically assume 
that temperature changes only in the direction perpendicu¬ 
lar to the surface (i.e. a plane-parallel approximation). This 
assumption allows a solution of the one-dimensional heat 
diffusion equation for each surface facet independently. It 
is well justified as long as surface features are significantly 
larger then the diurnal thermal skin depth, which varies from 
mm to dm ( Vokrouhlicky and Broz|1999 l. This assumption 
is no longer applicable for a high-resolution shape model 
and a full three-dimensional solution of the heat diffusion 
equation is required. Second, no shape is described to the 
required level of detail. So far, the best shape model is that 
of the asteroid Itokawa. The model in the best available res¬ 
olution consists of over 3 million facets, which corresponds 
to meter-sized surface features ( Gaskell et al.|[2C)06 |. 

As Golubov and Krugly (20121 pointed out, surface fea¬ 
tures of sizes comparable to the thermal skin depth could 
potentially have significant influence on the total YORP ef¬ 
fect. They considered a stone wall (an idealized boulder) lo¬ 
cated on the equator of a spherical asteroid and aligned with 
a local meridian. The wall was assumed to be high enough 
so that the heat would be mostly conducted in a transverse 
direction and the heat diffusion equation can be solved us¬ 
ing the one-dimensional approximation. They demonstrated 
that the emission from the surface of the wall can create 
a torque that will not vanish after averaging over the ro¬ 
tational period. Assuming a large number of such “walls” 
placed along the equator, the corresponding torque may 
be comparable to the torque arising from the global-shape 
asymmetry. 


Golubov et al. (2014 1 generalized the previous model by 


assuming spherical boulders. They studied a dependence of 
the torque on a number of free parameters of the problem. 


KdnU + eau'^ = 



d„u = 0 


Figure 1. The domain f! and boundary conditions of our prob¬ 
lem. The boulder is located on the top. The strips indicate the 
boundary r 2 , where the temperature is kept by a Dirichlet bound¬ 
ary condition. The temperature distribution inside the domain 
and on the surface Ti is then computed numerically as a solution 
of Eq. § with boundary conditions, Eqs. § to §. 


As expected, the resulting torque is lower than in a simple 
one-dimensional model; nevertheless, boulders can still af¬ 
fect rotational dynamics significantly. Authors assumed even 
mutual shadowing and heating of the boulders. 

The goal of this paper is to solve the heat diffusion 
equation in a realistic boulder. The problem requires a nu¬ 
merical solution in a general three-dimensional domain. We 
derive a formulation of the numerical problem in Section 
We discuss the magnitude of the torque induced by a single 
boulder in Section[^ We estimate the total torque that boul¬ 
ders contribute to the YORP effect on the asteroid Itokawa 
in Section]^ Finally, the results of our model and the impli¬ 
cations are summarized in Section]^ 


2 THE HEAT DIEEUSION EQUATION AND 
A WEAK FORMULATION OF THE 
PROBLEM 

Our problem may be specified as follows. We search for 
a temperature u(r, t) inside the boulder and its surround¬ 
ings, i.e. an unknown scalar function on a domain fl. The dif¬ 
ferential operator corresponding to the heat diffusion equa¬ 
tion (HDE) is: 

L = pCdt - V • AV , (1) 

where K denotes the thermal conductivity p the density, 
C the specific heat capacity of the material. The function u 
thus has to fulfill the relation: 

£(u)=0. (2) 

At the same time, we require the boundary conditions to be 
met at the boundary dQ. of the domain. 


2.1 Boundary conditions 

The boundary consists of two parts denoted Fi and r 2 , as 
shown in Figure[^ The boundary Fi represents the surface of 
the asteroid, the boundary condition is essentially an energy 
balance equation: 

KdnU + eau* = T, (3) 
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The thermal emission from boulders 3 


Notation 

Meaning 

u 

temperature 

u 

numerical solution 

c 

differential operator corresponding to 
the heat diffusion equation 

Nj 

basis functions 

W, 

weighting functions 

Uj 

coefficients of the decomposition of u 

K 

thermal conductivity 

P 

density 

C 

specific heat capacity 

€ 

infrared emissivity 

A 

hemispherical albedo 

(7 

Stefan-Boltzmann constant 

U) 

angular frequency 

p 

rotational period 

c 

speed of light 

<s> 

solar flux at the distance of the asteroid 

T 

flux absorbed by the surface 

L 

diurnal thermal skin depth 

e 

thermal parameter 


subsolar temperature 

n 

dimensionless pressure 

V 

latitudinal dependence of (fl) 

n 

computational domain 

Fi, Fa 

boundaries of the domain 

V 

position vector 

t 

time 

fi 

outward normal vector 

s 

body-Sun direction 

dn 

derivative along the (outward) normal 

P 

shadowing function 

V 

visibility function 

H 

Heaviside step function 

f 

force 

T 

torque 

e 

rotational axis 


asteroidal latitude 

1 

characteristic size of the boulder 

I 

moment of inertia of Itokawa 

N{i) 

differential size distribution of boulders 

7 

exponent of the size distribution 


Table 1. Notation used in the paper. 


where denotes a derivative along the normal, e the in¬ 
frared emissivity, a the Stefan-Boltzmann constant, the 
incoming flux absorbed by the surface. The boundary Fi is 
a non-convex surface, thus the radiative heat exchange also 
contributes to the total flux T (i.e. the self-heating effect). 
We denote the solar flux contribution Tq and the contri¬ 
bution of flux coming from visible parts of the surface — 
either the thermally emitted flux or the scattered flux — as 
.Fth and Tsc respectively. The total absorbed flux is the sum 
T = Tq -I- J^th + Tsc. These terms can be expressed as: 


T(^ = (1 — ■ n , 


= (1 — 


A) f dc 

Jri 


,4 COS a cos a 

'J, —7^- 

7r(r - fy 


udF'. 


(4) 

(5) 


-r- f cos a cos a 

Fac = (1- A) A^fi s - n ——=—j/dF , (6) 

Jr 'n(r-r'Y 


where A is the hemispherical albedo, >1> the flux of solar 
radiation, s the body-Sun direction, n the outward normal 


to the surface, r the position vector, a the angle between 
the local normal and the direction vector connecting points 
r and Y, jj, the shadowing function, v the visibility function 
(see Table . The prime denotes a value of a quantity at 
the point of surface element dF'. We assume the Lambert’s 
cosine law for the intensity of thermal emission and scattered 
radiation, hence the cos a' in equations ([^ and Q. 

The shadowing function p is defined on the surface Fi. 
The value of /i(r) equals 1 if the point r is insolated, 0 if it lies 
in the shadow. The visibility function v is defined on Fi x Fi 
space. We assign the value of function v{r, Y) to 1 if points f 
and r* have a visual contact, 0 otherwise. In most cases, the 
value of V is simply v(r, Y) = H(n ■ (r — r'))H{n' ■ {Y — r)), 
where H is the Fleaviside step function and n, n' denote the 
local normal at points rand d. 

The discretized boundary Fi is defined by a set Si of 
triangular facets. Integrals in equations and § can be 
therefore computed as sums, dF' —>■ Si. The values of 

the shadowing function fi and the visibility function u always 
correspond to whole facets. This restriction gives rise to an 
error; however, it can be estimated and limited substantially 
by choosing a high-resolution surface mesh. 

We also need to specify conditions on the boundary 
F 2 , which goes through the interior of the asteroid, closing 
the boundary dQ. It can be selected arbitrarily; we choose 
a shape corresponding to five walls of a block, which is a con¬ 
venient choice as we can simply set a zero-flux boundary 
condition: 


Kd„u = 0 . (7) 

This condition will be met as long as dimensions of the do¬ 
main fl are significantly greater than dimensions of the boul¬ 
der. The influence of the boulder can be considered negligible 
at large distances. At sides of the domain, the temperature 
will only change in the direction perpendicular to the sur¬ 
face, the dot product of the normal vector and the temper¬ 
ature gradient will therefore be null. At great depths, the 
temperature will be effectively constant, which means the 
temperature gradient at the bottom of the domain will be 
null, satisfying the boundary condition 0. 

We also need to specify an initial condition as the HDE 
is an evolution equation. However, we seek for a stationary 
solution that does not depend on a chosen initial condition. 
The choice of a condition will affect the speed of convergence 
only. The best available estimate of the solution is the lin¬ 
earized analytical solution in a half-space domain, derived 
in Appendix]^ hence we set: 

U = Utheoryi t = 0 . (8) 


2.2 A discretization 


We are going to solve the HDE (Eq. numerically, using 
a hnite-element discretization in space. In this approach, the 
function u is approximated by ( Langtangen|20d3 1: 


= u = ^UjNj 


(9) 


where Nj denote prescribed basis functions, Uj unknown co¬ 
efficients we search for and M corresponds to the number of 
vertices defined on the domain. Since u is only an approx¬ 
imation of u, applying the operator would generally yield 
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a non-zero result; 


( 10 ) 

nevertheless, we require the integral of all residua over the 
domain to be zero: 

[ C{u)W^An = Q, ( 11 ) 

Jn 

where Wi are suitable weighting (test) functions. This is 
called a weak formulation of the problem. In the Galerkin 
method, the test functions are simply the basis functions, 
Wi = Ni, so that: 


I £{u)Nidn = 0 . 
Jn 


( 12 ) 


express the corresponding matrices, because this is done au¬ 
tomatically by the FreeFem-l—I- code | Hecht|[2012 I. We use 
a conjugate gradient method for the matrix inversion, which 
is suitable for sparse linear systems. 

After a careful testing of our numerical method (see 
Appendix we choose the time step At = 10 ^P, where 
P is the period present in the insolation function P. The 
spatial step is controlled by the maximum volume Afl — 
10“® m^ of tetrahedra generated by the TetGen code (|Si| 


20061 . 


3 THE MEAN TORQUE CAUSED BY AN 
IRREGULAR BOULDER 


Essentially, this constitutes a system of M equations for Uj 
coefficients. 

In our case of the HDE (Eq. [^: 

[ pCdtuN.dQ. - / V • (A'Vu)Wdn = 0 . (13) 

Jn Jn 


The second term may be rewritten according to the Green 
lemma as: 


[ V ■ (K\7u)NidD. = - [ KVu-VNidn + i Kd^uNidF, 
Jn Jn Jan 

(14) 

which enables to incorporate the boundary condition easily, 
because we can express the normal derivative from boundary 
conditions ([^ and Q, that is KdnU = —eaiP + J- on Fi, 
KdnU = 0 on r 2 . 

Regarding the temporal derivative, we use a finite- 
difference discretization: 


dtu ~ 


At 


(15) 


and an implicit Euler scheme, so that we plug u" in the re¬ 
maining terms, whenever possible. The only exception is the 
non-linear radiative term, where we perform a linearization: 


/^n,m\4 / ^ n.m—1 \ 3 ^ n,m 

(ti ’ )U 


(16) 


and we employ an iterative method to find a solution; at 
the given time-step (denoted by superscript n) we find a 
sequence of solutions u™’’" to a linear problem, until — 
jg sufficiently small. The initial value u”’’^ can be 
selected arbitrarily, a common choice is 

In some cases, the iterative method does not converge 
and we thus introduce a relaxation parameter (. In each 
iteration, we find a preliminary solution m*’"* of the linear 
problem and we then assign a new value of u"’™ by taking 
a linear combination of the current and previous solutions: 

u = ;u*’ -I- (1 - Qu . (17) 


We achieved a convergence for all considered values of pa¬ 
rameters by selecting ^ = 0.6. 

The final equation is thus: 

[ / f^u^~^Nidn+[ K\/u^’”"-\/Nidn + 

Jn Jn Jn 

-b [ ecr(M"’’"“^)^M"’’"Wdr- [ Pdr = 0. (18) 

Jri Jri 

We actually need not to substitute for u from Eq. ([^ or 


The magnitude of a recoil force varies during a rotational pe¬ 
riod and a revolution of an asteroid around the Sun. A long¬ 
term effect of the force is therefore given by its time-averaged 
value. We follow the assumption of [Golubov and Kruglyi 
120121 and consider an asteroid on a circular orbit with zero 


obliquity. In fact, the eccentricity of Itokawa is e = 0.28; we 
address this issue in Section[3.5| Although the YORP effect 


depends on the obliquity in a non-trivial way (Capek and 
I Vokrouhlicky|[^04[ ) , the zero obliquity allows us to average 
the recoil force over a rotational period only. 

We consider the thermal emission and the scattered ra¬ 
diation. The direct radiation pressure has a negligible influ¬ 
ence ( Nesvorny and Vokrouhlicky||200^ K Again, we assume 
the Lambert’s cosine law for the intensity of scattered and 
emitted radiation. The recoil force from the surface element 
dS is then: 


dfth = - ‘^—u'^ndS, 

3 c 

d/ac = I^(s • n)ndS . 

3 c 


(19) 

( 20 ) 


The total torque caused by the boulder is given by the sur¬ 
face integral over the boulder: 


L. 


= / rxdf. 


( 21 ) 


The direction of the torque is generally different from 
the axis of rotation e. Both the direction and the magnitude 
of the torque depend on exact shape of the boulder. How¬ 
ever, even a symmetric boulder can induce a non-zero torque 
due to the lateral heat diffusion. The torque is caused by the 
asymmetry of emission from the eastern side and western 
side of the boulder. We anticipate the torque will therefore 
have a direction of the rotational axis e. 


3.1 The coordinate system and free parameters of 
the problem 


We choose a topocentric coordinate system centered on the 
studied boulder. The 2 axis has therefore a direction of a lo¬ 
cal normal, x axis is aligned with a meridian and y axis 
completes a right-handed orthogonal Cartesian system. 

We introduce quantities that help us to reduce a number 
of independent parameters of the problem. We define the 
subsolar temperature: 


M* 




( 22 ) 
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Figure 2. An example of the boulder shape which was used to 
compute the surface temperature distribution. 


the diurnal thermal skin depth: 


L = 


2K 

UJpC 


(23) 


where ui is the angular frequency of the asteroid, and the 
thermal parameter: 




(24) 


4\/27r 4 ecruj 

Numerical factors in these definitions arise from the deriva¬ 
tion of the analytical solution (see Appendix A|); however, 
some authors do not use them (|Lagerros|1996 Golubov and 
|Krugly|201^ . 

If we neglect self-heating terms, the heat diffusion equa¬ 
tion § and its boundary condition § can be rewritten 
using dimensionless variables ^ = r/L, 95 = uit, v = u/m* as: 


1 . dv 

47 r On • V V =s-n, 


(25) 

(26) 


where Vj, Aj is the gradient and the Laplacian with re¬ 
spect to the variable The only independent parameter in 
these equations is the thermal parameter 0. However, the 
boundary condition must hold for all £ dQ. If £ is the 
characteristic size of the boulder, then the problem of finding 
a dimensionless temperature v has actually two independent 
parameters — the thermal parameter 0 and the dimension¬ 
less size £/!/. In the following, we select the characteristic size 
as the square root of the base area of the boulder, £ = y/S. 


3.2 The dimensionless pressure 


In our model, the shape of the boulder can be arbitrary. 
As a special case, we can choose a high wall and compare 


our results to the model of Golubov and Krugly (20121. We 


started with such an idealized “boulder” to determine the 
influence of the self-heating effect and the absorption of ra¬ 
diation on the recoil force (see Appendix E- Nevertheless, 
we then selected a boulder of realistic irregular shape for 
our computations, as shown in Figure The shape was ob¬ 
tained by a 3D scanning of a randomly selected boulder. In 
this case, we do not take into account the self-heating ef¬ 
fect nor the influence of absorption though; this decision is 
justihed in Appendix [B| as well. 

We considered different values of the thermal conduc¬ 
tivity for the studied boulder and for the surrounding layer 
of regolith, as demonstrated in Figurej^ Thermal properties 



Figure 3. A sectional view of the 3D mesh, created by the Tet- 
gen code. Thermal parameters of the boulder are varied in our 
model. Nonetheless, we assume the boulder has generally differ¬ 
ent thermal conductivity than the surrounding layer of regolith. 
Therm al properties of the regolith as given by [Farinella et ah] 
( |l998| are as follows: the density p = 1500 kg m“^, the thermal 
conductivity K = 0.0015 Wm“^ and the specific heat ca¬ 
pacity C = 680Jkg“^K“^. We chose the hemispherical albedo 
A = 0.1, the infrared emissivity e = 0.9, and the rotational period 
equal to P = 12.1 hours. 


of the regolith are taken from [Farinella et al.| ( [T9M| , prop¬ 
erties of the boulder are determined by values of thermal 
parameter 0 and the skin depth L. 

We should stress the importance of the non-linearity of 
the problem. We derived a linearized analytical solution of 
the heat diffusion equation in a half-space domain (see Ap¬ 
pendix 0 - where we deal with the non-linear term by 
substituting Mq -b 4mq5u, where uo is a constant, 5u is the 
change of temperature. The same term appears in the ex¬ 
pression for the recoil force ( |19| ) from a thermal emission. 
In case of a symmetric boulder, the complete linearization 
would lead to identically zero mean torque. Therefore, we 
solve a non-linear problem by an iterative method, as de¬ 
scribed in Section!^ 

The solution of the HDE is a time-dependent tempera¬ 
ture distribution in the boulder, particularly on its surface. 
It follows that we can determine the recoil force the boulder 
exerts; the force is given by the formula (19l. However, it is 
convenient to introduce the dimensionless pressure: 


H = 


2 1 
35 


L 


, riy dr , 


(27) 


where Uy is the y-th component of the local normal, S is 
the base area of the boulder. The dimensionless pressure 
allows us to compare the magnitude of the tangential force 
for different sizes of the boulder. 

The projection of the total torque to the rotational axis 
is then given by: 


-HSrx, 


(28) 


where ry is the distance of the boulder from the rotational 


If we consider a wall aligned with a local meridian, 
which face of area S has a constant temperature u, our def¬ 
inition (271 is then reduced to: 

which is equivalent to the definition of a dimensionless pres¬ 
sure by Golubov and Krugly (20121. Our definition is there- 
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l/L 



l/L 


Figure 4. Computed values of the mean dimensionless pres¬ 
sure (n) as a function of the dimensionless boulder size ijL for 
two different values of the thermal parameter © = 5.38 and 1.70 
(notice that vertical axes have different scalings). The dark curves 
correspond to the boulder rotated by 0°, 90°, 180° and 270° 
around the vertical axis. We notice that all curves approach a zero 
for IjL 0 and they exhibit a maximum for £ ~ L. For © = 1.70 
there seems to be an inflection at about £ ~ 0.1 L, which is absent 
for © = 5.38. The bright curve is the arithmetic mean of the dark 
curves. We assumed a simpler “shadowing” model in this case 
(but cf. Figure [bT] |. 


fore analogous, with the exception of the area S being the 
projection to the plane xy rather than the plane xz. 

The dimensionless pressure varies during a rotation. We 
can obtain a measure of the long-term effect by averaging 
over one rotational period; to this point we introduce the 
mean dimensionless pressure: 

{n) = ^iydt. (30) 


3.3 The total pressure exerted by a set of 
variously oriented boulders 

The mean dimensionless pressure (II) as a function of the di¬ 
mensionless size IjL varies significantly for different shapes 
of the boulder, or even for different orientations of the same 


bonlder. It is evident that the limit of very high conductivity 
(that is f/L —>■ 0) leads to a zero dimensionless pressure II 
for all shapes of a boulder. In snch a case, the boulder is 
isothermal and therefore emits the same flnx to the western 
and eastern directions, resulting in a null torque. 

The limit of the mean dimensionless pressure for zero 
thermal conductivity {l/L —>■ oo) differs from boulder to 
boulder. The conductive term in the energy balance equa¬ 
tion © is negligible and the temperature at a given point of 
surface is determined by the immediate balance between in¬ 
coming and outgoing radiant flux. Since we solve the HDE 
in a single boulder only, we need to obtain a torque (as 
a function of a boulder size) that would represent a set of 
all boulders on the surface. It is reasonable to assume that 
the boulders are randomly oriented on the surface. Although 
some orientations of boulders seems to be preferred on cer¬ 
tain parts of the surface of Itokawa ( Miyamoto et al.|2007 l, 
we anticipate that no orientation prevails on a global scale. 
The total torque induced by boulders will therefore vanish 
in the limit of zero conductivity. For that reason, we demand 
the mean torque (II) to approach zero as well. However, in 
the general case of an asymmetric boulder, the mean dimen¬ 
sionless pressure will approach a non-zero value. 

We have several options how to resolve this issue. For 
instance, we can restrict ourselves to symmetric boulders 
only. If the boulder is symmetric with respect to the plane 
of the local meridian, the mean pressure will vanish in the 
limit case. Nonetheless, we want to keep the universality of 
our model and use irregular asymmetric boulders. In this 
case, we can compute the mean pressure (H) for several ori¬ 
entations of the boulder and then take the average of these 
values. Another possibility is to calculate the mean pressure 
for a single orientation and subtract the pressure in the limit 
of zero conductivity. We employed the former option. 

Figurej^shows the mean dimensionless pressure for sev¬ 
eral orientations of the studied boulder and the averaged 
values. We assumed the boulder lies on the equator of an 
asteroid. We see that the averaged values approach zero in 
the limit of zero conductivity, as expected. 


3.4 A dependence of the mean pressure on 
asteroidal latitude 

For an asteroid with zero obliquity, the body-Sun vector s 
has Cartesian coordinates: 

s = (sin cos HA, — sin HA, — cos cos HA), (31) 

where HA is the hour angle and is the asteroidal latitude 
(defined as sini? = e • n). The dependence of the dimension¬ 
less pressure H on the hour angle HA vanishes after averag¬ 
ing over a rotational period, the dependence on i? remains. 

Assuming we can separate variables £, we can decom¬ 
pose the mean pressure (H) as: 

{n)(f,r?) = P(f)V(r?), (32) 

where P{£) = {n)(f, 0). The function V(d) constitutes a lat¬ 
itudinal dependence and is normalized such that V(0°) = 1. 
It obviously depends on the shape of a boulder. For this 
test, we chose an approximately hemispherical boulder, be¬ 
cause it is axially symmetric and therefore one latitude is 
not preferred over other values. 

We show the computed values of the function V(i?) in 
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Figure 5. The dependence of the mean dimensionless pres¬ 
sure (n) on the asteroidal latitude ^ of the boulder. 


Figurej^ It can be approximated by a function cos at), where 
a = 0.653±0.004 is a parameter determined by a least square 
fit. The mean pressure (II) is maximal for the boulders lo¬ 
cated on the equator and does not drop below 50% of the 
maximum for t? = 80° latitude. The known dependency of 
the mean pressure on the latitude allows us to compute the 
torque induced by a boulder on any given point of the sur¬ 
face, which we shall use in the next section. 


3.5 An influence of the elliptical trajectory 

So far we assumed the asteroid orbits on a circular trajec¬ 
tory. This assumption allowed us to ignore any seasonal ef¬ 
fects and average the torque over a rotational period only. 
However, the eccentricity of Itokawa is e = 0.28; it is not 
clear whether the ellipticity of the trajectory affects our 
model significantly. We thus compute the mean dimension¬ 
less pressure (H) in several “discrete” points on the orbit 
separated by a time step equal to 1/30 of the orbital period. 
We assume a high-conductivity case here. Nevertheless, the 
thermal parameter 0 will vary over time; it depends on the 
distance from the Sun as 0 ~ r5, so it will reach 0 = 3.29 
in the perihelion and 0 = 7.79 in the aphelion. We restrict 
ourselves to a single value of the dimensionless size i/L = 1 
and we average the result over four different orientations of 
the boulder, as before. Then, we obtain the time-averaged 
value by simply taking an arithmetic mean of the computed 
values. 

This approximate average over elliptical orbit is to be 
compared with the mean dimensionless pressure correspond¬ 
ing to the circular orbit, which we obtained already in Sec¬ 
tion |3.3[ The result can be seen in Figure We see that 
even for the considerable eccentricity of Itokawa, the differ¬ 
ence seems negligible. The value differs from our previous 
result by less than 5 %. Therefore, the approximation of the 
circular trajectory is well justified. 


Figure 6. The values of the mean dimensionless pressure (11} 
from the perihelion (t = 0) to aphelion (t = Pl2). The values are 
computed for the high-conductivity case (0 = 5.38 at the distance 
of the semimajor axis) and for I = L, averaged over four different 
orientations of the boulder. One can readily see that the difference 
between the time-averaged value and the value of corresponding 
circular orbit is negligible, given other uncertainties of our model. 


4 THE ANGULAR ACCELERATION OF 
ASTEROID (25143) ITOKAWA 

In the following, we focus our attention on the asteroid 
(25143) Itokawa. First, we estimate the total number of 
boulders on the surface and then compute how the ther¬ 
mal emission from boulders alters the angular acceleration 
predicted by global-shape models of the YORP effect. 

Existing models of the YORP effect usually assume 
a normal direction of the recoil force. However, for non- 
convex asteroids the force can be influenced by the absorp¬ 


tion of radiation emitted by the surface (Statler 20091. In 


the previous section, we demonstrated that a surface feature 
can alter the recoil force as well. We pointed out that the lat¬ 
eral heat diffusion through boulders leads to an emergence 
of the tangential component of the recoil force. The pres¬ 
ence of surface features also changes the normal component. 
While a complete solution would require solving the heat 
diffusion equation in the whole asteroid, including boulders, 
we neglect the change in the normal component and we solve 
for the tangential component separately. 

The torque generated by a single boulder was discussed 
in previous chapter We now place a large number of such 
boulders on the shape model of Itokawa and calculate the 
torque they induce. The total YORP torque and correspond¬ 
ing change of the angular velocity of the asteroid is then ob¬ 
tained by adding our result to the result of the global-shape 
model of the YORP effect. 


4.1 The torque induced by boulders 

We demonstrated the emergence of the tangential compo¬ 
nent of the force, which is of our interest. Therefore, we 
consider a direction of the force perpendicular to the lo¬ 
cal normal n, regardless of the location on the surface of 
the asteroid. Although the direction depends on the shape, 
we discussed that the overall effect corresponds to the force 
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perpendicular to the rotational axis e. In follows that the 
direction of force is e x n. 

The magnitude of e x n is proportional to cosi9. How¬ 
ever, we studied the dependence of the torque on asteroidal 
latitude in Section [3.4[ We showed that the mean torque as 
a function of latitude is approximately cos0.653i?. The re¬ 
coil force induced by single boulder is therefore proportional 
to the term: 




cos0.653i? 

cos-!? 


e X n. 


(33) 


Our goal is to compute the total torque r induced by the 
boulders; to be more precise, its projection r • e to the rota¬ 
tional axis. Let N{l)Al be the number of boulders in the size 
interval [1,1 -\- Al). We can write the size distribution N{1) 
as: 


N{t:)Al = Noq{1)AI, 


(34) 


where Nq is the total number of boulders on the surface and 
q{1)AI is the probability that a randomly selected boulder 
has a size in the interval -|- df?). We compute the to¬ 
tal torque induced by boulder as follows. We select a facet 
of Itokawa shape model randomly. Assuming boulders cover 
Itokawa uniformly, the probability of selecting each facet is 
proportional to its area. Now, we generate a random size I 
with the probability distribution q, using the inverse sam¬ 
pling method. The torque induced by the boulder on the 
selected facet is added to the total torque. The final num¬ 
ber of generated boulders is given by the value No from 


Eq. (34l. Because the number of boulders is very high, the 


total torque basically does not depend on the realization of 
boulder distribution on the surface. 

It is clear that the boulder size distribution N[t)Al con¬ 
stitutes an important parameter for a calculation of the total 
torque. We thus discuss the size distribution of boulders on 
Itokawa in the following section. 



Figure 7. The image ST_2563607030-v jSaito et al.||2010[ | with 
highlighted boulders from which we derived their size distribution, 
used for the computation of the total torque. 



I [m] 


4.2 The observed size distribution of small 
boulders 

In order to obtain the torque caused by boulders, it is nec¬ 
essary to find out the total number of boulders and their 
size distribution. The differential size distribution of boul¬ 
ders larger than 5 m on the whole surface of Itokawa can be 
approximated by a power law ( Saito et al.|2006 1 : 

N{£)Ae ^ 1.3 X 10® [f]“® ®df. (35) 

Surface images taken by the Hayabusa spacecraft revealed 
that this power law can be extrapolated down to sizes of 
0.1 m on certain parts of the surface ( Miyamoto et al.|2007 l, 
although the slope of a log-log graph falls off significantly 
for smaller sizes. However, other parts of the surface clearly 
show a different topography. Furthermore, an extrapolation 
of the above-mentioned size distribution down to 1 mm is 
clearly unacceptable; boulders of sizes between 1 mm and 
0.1 m alone would cover about 4 x 10^ nP of the surf ace area, 
but the surface of Itokawa is only 3.93 x 10® m^ (Demura 
et al.|2006 l. 


Therefore, we sought for a different size distribution 
of small pebbles. We estimated their size distribution from 
close-up images taken by the Hayabusa during its descend, 
namely the images ST_2563537820_v and ST_2563607030_v 


Figure 8. The histogram of small boulder sizes on the sur¬ 
face of Itokawa, constructed from images ST_2563537820_v and 
ST_2563607030_v. The three power-law fits correspond to the av¬ 
erage fit and left/right l-tr uncertainty, determined by a Monte- 
Carlo method. Two rightmost bins lie off the linear part of the his¬ 
togram and were not included in the least-square method dataset. 
The values of the multiplicative constant and the exponent of the 
power law are shown in Eq. (|36[|. 


from Saito et al. (20101 dataset. The resolution of these im¬ 


ages is 7mm/pixel and 6mm/pixel, respectively (Miyamoto 
et al. ][2007| ), which allows us to find distinct boulders only 
few centimeters in size. Identified boulders are shown in Fig¬ 
ure We constructed a histogram of sizes (see Figure 
from which we inferred the size distribution: 


N{£) A£ = (14 ± 9) X 10® d^. 


(36) 


We assume this power law can be extrapolated to millimeter¬ 
sized pebbles. Considering the uncertainties, the area cov¬ 
ered by boulders of sizes between 1 mm and 1 m would be 
between 19 % and 32 % of the total surface area. This finding 
seems to be viable. 
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4.3 A comparison of the YORP torque by 
boulders and by the global shape 

The global-shape YORP effect model of the asteroid Itokawa 
predicts a significant rotational deceleration, which is incon¬ 
sistent with the observed acceleration. As mentioned above, 
the lateral heat diffusion through boulders induces an addi¬ 
tional torque, which affects the angular velocity. 

Let the magnitude of the total torque generated by 
boulders be r. The asteroid will undergo the rotational ac¬ 
celeration: 


duj 

dt 


dcu\ 
dt ) 


+ 


global 


I ’ 


(37) 


where (da;/dt)giobai is the prediction of the global-shape 
YORP model, I = 7.77 x 10^"* kg nd is the moment of inertia 
of Itokawa ( Scheeres et al.|[2007 1. 

The global-shape model of the YORP effect predicts 
a rotational deceleration —(2 to 6) x 10“^radday“^, de¬ 


pending on the resolution of the shape model | Breiter et al. 


2009). In order to determine the torque induced by boul¬ 


ders, it is necessary to select values of the parameters — 
namely the thermal parameter 0 and the thermal skin 
depth L. We adopted following material properties: K = 
2.65 Wm-i K-i, C = 680 Jkg-^ K-\ p = 2700kgm-^. To¬ 
gether with orbital parameters of Itokawa, this yields the 
thermal parameter 0 = 5.38 and the thermal skin depth 
L — 0.141m. Utilizing the size distribution of boulders de¬ 
rived in Section 14.21 we obtain the result: 


(38) 


= (1.20 ± 0.11) X 10“’’radday“^ . 


I le=5.38 

As the thermal conductivity seems to be the most uncer¬ 
tain parameter, we also computed the torque for a lower 
value, K = 0.26Wm“^K“'^ (keeping other parameters un¬ 
changed). In such a case, the thermal parameter is 0 = 1.70 
and the thermal skin depth L = 0.0446 m. The correspond¬ 
ing torque is then: 


= (4.8 ± 1.2) X 10"’’radday“^ . 


(39) 


/ le=i.7o 

The probability distribution for both cases is shown in Fig¬ 
ure |§1 

For the sake of comparison, we can refer to the result of 


global-shape models of|Lowry et al. 

(12014), (-1.80±1.96)x 

10 ^ rad day , or 

Breiter et al.| 

2009), -(2.5 to 5.5) x 


10“ rad day” . We notice that this torque is of the same 
order as our result, but has an opposite sign. The torque in¬ 
duced by boulders and the torque from the global asymme¬ 
try could effectively cancel out, resulting in the change of an¬ 
gular velocity much smaller than predicted by global-shape 
models. As the observed ang ular acceleration o f Itokawa is 
(0.35±0.04) X 10”^radday”^ (Lowry et al. 2014), our model 


presents an alternative explanation of the observed acceler¬ 
ation without any need for a non-uniform density distribu¬ 
tion. 


5 CONCLUSIONS 

In this paper, we presented a detailed numerical model of 
the local YORP effect induced by a boulder or a set of boul¬ 
ders. The three-dimensional heat diffusion equation in the 
boulder was solved using the finite element method. Unlike 




3 4 5 6 7 8 9 

t/I [10“^ rad day“^] 


Figure 9. The contribution of boulders to the total angular accel¬ 
eration of Itokawa, computed for two different values of the ther¬ 
mal parameter 0 = 5.38 and 1.70 and the corresponding thermal 
skin depth L. Each histogram was constructed by a Monte-Carlo 
method; individual runs correspond to different power laws of 
the boulder size distributions. Note that horizontal axes differ in 
each figure. Eqs. | |38^ and | |39[ l show the average value and l-tr 
uncertainties. 


the finite difference method, the finite elements have basi¬ 
cally no restriction on the shape of a domain, allowing us 
to solve the heat diffusion equation in the boulder of a re¬ 
alistic shape. Furthermore, we assumed the studied boulder 
has a different thermal conductivity than the surrounding 
regolith layer. 

Our boulder had a general asymmetric shape, so it ex¬ 
hibited a non-zero torque even in the limit of zero thermal 
conductivity. However, this torque depends on the orienta¬ 
tion of the boulder. In order to obtain an average torque 
representing a set of randomly-oriented boulders, we com¬ 
puted torques for several orientations and then the average 
of these values. We verified that the averaged torque ap¬ 
proaches zero in the zero conductivity limit. 

The non-zero torque arises from the asymmetry of the 
emission (averaged over the rotational period). There are 
two reasons for the emission asymmetry. The first is the 
asymmetry of the boulder shape. Indeed, Rubincam (20001 
demonstrated an emergence of the YORP effect on a toy 
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model of a spherical asteroid with two wedges attached to 
its equator. The torque created by the emission from the 
vertical side is greater in magnitude than the torque cre¬ 
ated by the emission from the inclined side, thus resulting 
in a non-zero total torque. 

The second reason for the emission asymmetry is the 
lateral heat diffusion through the boulder (as in [Golubov 


and Krugly||2012 1 . We can imagine a boulder on the equa¬ 


tor. In the morning, the eastern side of the boulder is heated 
up and the boulder exerts a recoil force of western direction. 
Provided the width of the boulder is comparable to the ther¬ 
mal skin depth, the heat diffusion contributes to heating of 
the western side in the afternoon. The emission from the 
western side is therefore more intense. The recoil force has 
an eastern direction and exceeds the force from the eastern 
side in magnitude, thus creating a non-zero mean force of 
eastern direction. The corresponding torque causes an an¬ 
gular acceleration of the asteroid. 

The global contribution of the shape asymmetry of 
boulders to the YORP effect is likely to be null, because of 
the very large number of boulders on the surface. In contrast, 
the lateral heat diffusion leads to the torque with a direc¬ 
tion of the rotational axis, thus accumulating over individ¬ 
ual boulders. Even though the torque generated by a single 
boulder is tiny, the overall effect can be comparable to the 
global-shape YORP effect, if there is a sufficient amount of 
boulders, of course. 

We showed the maximum pressure is exerted by boul¬ 
ders of sizes comparable to the diurnal thermal skin depth L. 
Figure 1^ shows the dependence of the pressure on the size 
of the boulder. However, this graph does not reflect the ac¬ 
tual contribution of boulders of different sizes to the total 
torque. The decisive factor is the exponent 7 of the power- 
law size distribution. If 7 > 2, smaller boulders will exert a 
higher torque (compared to the Fig.|^. In our case, the over¬ 
all torque is generated mostly by boulders of sizes between 
0.11/ and L. 

The general approach allowed us to compare our three- 
dimensional model with the one-dimensional model of IGol-l 


ubov and Krugly (20121. In case of a symmetric boulder, we 
confirmed that the torque vanishes in the limits of high con¬ 
ductivity and zero conductivity. We showed that the maxi¬ 
mum torque appears for t L. 


Unlike Golubov and Krugly (20121, we found positive 
values of the torque for all parameter ranges in the case of 
a symmetric boulder. An asymmetric boulder could produce 
a negative torque, but after averaging over orientations the 
resulting torque is again strictly positive. Even jGolubov and) 
Krugly (20121 realized the torque is mostly positive, and 


proposed a possibility of an equilibrium between the global- 
shape torque and the torque induced by boulders, resulting 
in a null total torque; they suggested this could be the case of 


the asteroid (25143) Itokawa. However, Lowry et al. (20141 
detected a positive change in angular velocity of Itokawa, 
which means this asteroid is not in such equilibrium state. 

Of course, our model contains a number of free param¬ 
eters that can change the magnitude of the torque signifi¬ 
cantly. The crucial factor is the total number of boulders and 
their size distribution. We realized that the size distribution 


of the surface of Itokawa and we extrapolated it for sizes of 
indiscernible pebbles. We also assumed that other parts of 
the surface have the same size distribution of boulders. If we 
neglect an interaction between boulders (mutual shadowing 
or thermal irradiation), which seems to be reasonable given 
the separations of boulders seen on Figure the torque is 
directly proportional to the number of boulders. 

The choice of the lower limit of the size distribution 
might be also a disputable parameter. Although the magni¬ 
tude of the torque induced by a boulder approaches zero as 
the size of the boulder approaches zero, even sub-millimeter 
pebbles could have a non-negligible influence on the total 
torque. However, it is doubtful whether so small particles 
can be considered as boulders, or whether they form a uni¬ 
form layer of matter. We selected the lower limit 1 mm. 

The shape of the studied boulder is another key factor 
of our model. We selected a boulder of a realistic irregular 
shape, but it was selected ad hoc. There are two particularly 
important properties of boulders: their height and flatness. 
The higher the boulder is, the greater torque it likely in¬ 
duces. If the sides of the boulder are perpendicular to the 
surface, the lever arm of the torque is maximal; the lower 
the slope of sides, the lower lever arm. 

We finally applied our model to the case of the aster¬ 
oid (25143) Itokawa. We showed that boulders could in¬ 
duce a torque that would cause the angular acceleration 
of the order 10“^ rad day We realized there was a sig¬ 
nificant uncertainty; nevertheless, we clearly demonstrated 
that the emission from boulders is capable of producing the 
torque comparable in magnitude to the global-shape YORP 
effect. Our result is consistent with the observed acceler¬ 


ation (Lowry et al. 20141 and presents an alternative and 


viable explanation of the discrepancy between the observed 
acceleration and the acceleration predicted by global-shape 
models. 


6 FUTURE WORK 

We postpone the following topics for future work. We as¬ 
sumed a rotational axis perpendicular to the orbital plane, 
therefore we need not consider the orbital movement and 
the torque is averaged over the rotational period only. Luck¬ 
ily, the obliquity of Itokawa is approximately 178° (Demura 


of large boulders on Itokawa of Saito et al. (2006 1 cannot be 
extrapolated to centimeter sizes. Having no better alterna¬ 
tive, we estimated the size distribution from close-up images 


et al.|[2006 l, which is very close to the perpendicular state. 
Considering a general direction of the rotational axis, the 
insolation will change during the revolution about the Sun 
and the seasonal changes of temperature will occur. We ex¬ 
pect a seasonal variant of the studied effect to appear on 
boulders whose size is comparable to the seasonal thermal 
skin depth. 

We already discussed that the diurnal torque has a di¬ 
rection of the rotational axis, as it is caused by the asymme¬ 
try of emission between the western and the eastern part of 
the boulder due to the lateral heat diffusion. Following the 
same principle, the seasonal torque could be caused by the 
asymmetry between the northern and the southern part of 
the boulder. The direction of this torque would therefore be 
perpendicular to the rotational axis. We thus anticipate the 
seasonal component of the effect will not affect the angular 
velocity of the asteroid, it will only cause an evolution of the 
obliquity. 
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In our model, we assumed a shadow is casted by the 
boulder; however, we did not take into account shadows 
casted by global-shape inconvexities of Itokawa. The same 
goes for the self-heating effect, which is also considered only 
locally. To resolve this issue, it would be necessary to de¬ 
termine the insolation function for each facet of the shape 
model separately and then solve many three-dimensional 
heat diffusion problems. 
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APPENDIX A: A LINEARIZED ANALYTICAL 
SOLUTION OF THE ONE-DIMENSIONAL 
HEAT DIFFUSION EQUATION 


The general three-dimensional heat diffusion equation with 
a non-linear boundary condition in an irregular domain has 
no analytical solution. To find the temperature distribution, 
we must employ a numerical approach such as the finite 
element method. Nevertheless, it is useful to derive an ana¬ 
lytical solution for a simplified case of a half-space domain, 
which allo ws to reduce t he problem to one spatial dimension 
only, as in Capek| ( |2007[ ) . The solution can be used as a test 
for our numerical model and also as a Dirichlet boundary 
condition (see Figure [^. 

Suppose the Sun illuminates an (infinite) plane z = 0 
and a half-space z > 0 represents the domain Q. We seek for 
the temperature u as a function of the depth z and time t, 
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solving the heat diffusion equation: 


ad^u — dtu = 0, (Al) 

with boundary conditions: 

—Kdzu(0, t) + eau^{0, t) = T{t) , (A2) 

9zu(oo,t)=0, (A3) 

where a = ^ is the thermal diffusivity, K the thermal con¬ 
ductivity, p the density, C the specific heat capacity, e the 
infrared emissivity, a the Stefan-Boltzmann constant and 
J-{t) the incoming radiant flux. The first boundary condi¬ 
tion is the energy balance equation and the second one is 
necessary for a uniqueness of the solution — it eliminates so¬ 
lutions where the temperature rises ad infinitum as 2 : —>■ oo. 

The radiant flux is a periodic function, therefore we can 
represent it as a Fourier series, J-{t) = X]^_oo We 

look for a stationary solution, which we can represen t by 
a sum u{z,t) = . Substituting into (All 

and applying the constraint (|A3l) we obtain the solution: 


uo(z) = ao , 
u„(z) = , 

u-„(z) = , 


(A4) 

(A5) 

(A6) 


where /3„ = \/ . Notice that /li is the reciprocal of the 


thermal skin depth L introduced in Section 3.1 hence the \/2 


factor. We determine constants a„ from the boundary con¬ 
dition (A21. Here we encounter problems with the non-linear 


term A substitution of the Fourier sum for the tempera¬ 
ture u would lead to a set of non-linear equations, which — 
unsurprisingly — does not have an analytical solution. Nev¬ 
ertheless, we can solve the problem analytically under the 
assumption that the changes of the temperature are signifi¬ 
cantly smaller than its absolute value. In such a case, we can 
linearize the fourth power ~ uq -|-4uo The 


boundary condition (A2 I then composes a set of linear and 


separated equations for the coefficients o„. We immediately 
see the solution for the constant term: 


Wo 


(A7) 


It is convenient to introduce auxiliary parameters that help 
us to get rid of complex numbers in the solution. First, we 
present the thermal parameter 0„ of the n-th mode: 

K/3n 


On = 


4ecrUQ 


(AS) 


Once again, 0i corresponds to the thermal parameter 0 
defined in Section [3.1| Second, let the phase shift ipn of the 
n-th mode be: 


tan tfin = — 


On 


On + 1 


sgnn. 


(A9) 


Since the solution is periodic, we can choose the initial time 
arbitrarily. Therefore, we choose the insolation function in 
the form of cosine series: J-{t) = (1 — A)4? H(costJt), where 
H(x) = X for X ^ 0, S(x) = 0 for 2 : < 0. First four Fourier 
modes of this function are: 

/II 2 2 

J^(t) ~ (1 —A)<E> ( - ^ cos Lot -\ - cos 2LOt - cos 4:LOt 

\n 2 3n IStt 

(AlO) 



d/L 

Figure Bl. The mean dimensionless pressure {11} as a func¬ 
tion of the dimensionless width d/L for four different variants of 
our model for an idealized boulder: i) shadowing only, ii) factor, 
iii) self-heating, and iv) a complete model. See text for a detailed 
explanation. 


Now, we can write the solution u of the problem in a 
simple form: 


u{z,t) 



Tn e cos (nuot — finZ + iPn) 

2efjWQ yJ20‘n + 20n T 1 

(All) 


We see that for each cosine term in the series of the inso¬ 
lation function F there exists a corresponding cosine term 
in the solution u, but with some phase shift. For the sur¬ 
face temperature in particular, the offset of the n-th Fourier 
mode is equal to ipn, defined above. 


APPENDIX B: VARIANTS OF THE 
INSOLATION FUNCTION AND 
A COMPARISON WITH THE EXISTING 
ONE-DIMENSIONAL MODEL 


The influence of topographic features on rotational dynam¬ 
ics of an asteroid has been studied by [Golubov and Kruglyi 
(20121. Their model uses a one-dimensional approximation. 


which implicitly corresponds to a “wall” of an infinite height. 
Such an infinite domain cannot be used in our numerical 
model; therefore, we used the wall the height of which is 
about three times the width. The vertical heat diffusion 
along the wall is therefore significantly reduced (in com¬ 
parison to the boulder studied in Section]^. Furthermore, 
we employed four variants of our model, called i) shadow¬ 
ing, ii) factor, iii) self-heating and iv) complete. They can 
be understood as sorted by a degree of their completeness. 

The first variant is called the “shadowing” model. It 
takes into account shadows casted by the boulder only. The 
self-heating effect is ignored, the influence of absorption on 
the direction of the recoil force is neglected as well. The 
insolation function is simply: 


= (1 — A)4?/rs • ft. 


(Bl) 


The second model, which we call the “factor” model, 
accounts for the self-heating by including simply the factor 2 
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to the solar flux. The corresponding insolation function is 
therefore: 


T = 2{l-A)^fj.s-n. (B2) 

This model has been used by [Golubov and Kmgly| ( |2012| ) . 

The third “self-heating” model computes the self¬ 
heating contribution directly by a calculation of the thermal 
emission and the scattered flux as described in Section 

T = {1 — A)^lj,s-n + Tth + , (B3) 

where J-tb, J'ac are the fluxes defined by Eqs. §. 

We called the fourth model “complete”, as it contains 
the direct computation of the self-heating, but also the in¬ 
fluence of the absorption on the direction and magnitude 
of the recoil force. We assume that the thermally emitted 
radiation falling on the surface is absorbed and does not 
contribute to the torque. The recoil force is then given by: 


dS" 


ir-L 


COS a cos a r — r 
7 r(f-f')2 ||T'-r|| 


vdV' ] (B4) 


The mean dimensionless pressure (II) (see Section 3.21 


as a function of the dimensionless width d/L of the wall is 
shown in Figure |BT| We notice certain common properties of 
the functions. All of them have a global maximum at d ~ L, 
but there is also either a secondary local maximum or an 
inflection at d ~ O.IL. These properties are evident for an 
irregular boulder as well, as we pointed out in Section [^ 
The maximal value of the “factor” curve is ~ 0.014, 
which is a result consistent with the findings of [Golubo^ 
and Krugly (20121. We also see that this model is a remark¬ 


ably viable approximation of the “self-heating” model. We 
thus confirm that the multiplication of the solar flux by the 
factor 2 leads to a similar outcome as the inclusion of the 
self-heating effect (which is much more computationally de¬ 
manding). 

Let us compare the “shadowing” and “self-heating” 
models. We see a notable property of the local YORP ef¬ 
fect: the self-heating effect causes an increase of the mean 
pressure by ~ 50 % in our model. The YORP effect induced 
by boulders therefore qualitatively differs from the global- 
shape YORP effect, where an inclusion of the self-heating 
effect leads to a decrease of the torque magnitude for most 
cases ( Rozitis and Green|2013 l. 

Finally, we compare the “shadowing” and “complete” 
model. We observe that the “complete” model is well ap¬ 
proximated by the “shadowing” model in the vicinity of the 
global maximum, although it differs significantly for smaller 
boulders. Nevertheless, this result allowed us to proceed with 
the simple “shadowing” model. Regarding the conclusions of 
this paper, a computation with the complete model would 
lead to very similar results. 


APPENDIX C: THE NUMERICAL 
UNCERTAINTY OF THE FINITE ELEMENT 
METHOD 

The discretization of both space and time inevitably intro¬ 
duces a numerical uncertainty to the solution. In order to 
estimate this uncertainty and find its upper bound, we uti¬ 
lize the linearized analytical solution of the HDE, derived 
in Appendix We select a simple block as a domain for 



Figure Cl. The maximal absolute difference of the numerical 
and analytical solution, plotted as a function of the discretization 
parameters — the maximum volume Afl of tetrahedra and the 
time step At. The limit S of the iterative procedure is constant, 
5 = 10-'‘. 


our numerical model, i.e. we use a flat surface without any 
boulders or other surface features. 

We are aware that the analytical solution is valid only 
if the amplitude of temperature is significantly smaller than 
its mean value. That essentially corresponds to high values 
of the thermal conductivity K (more precisely, of the ther¬ 
mal parameter 0). Namely, we select such parameters that 
the average temperature is ~ 180 K and the changes of tem¬ 
perature reach up to ~ 6K. In such a case, the linearized 
analytical solution is a sufficiently good approximation. 

The numerical uncertainty is influenced by the selection 
of three parameters: 


(i) The spatial discretization parameter Afl, which cor¬ 
responds to the upper bound for the tetrahedra volume in 
the three-dimensional mesh. We use this parameter during 
the generation of the mesh by the TetGen code | Si|2006 1 . 

(ii) The time step At. If P denotes the rotational period, 
then P/At time steps are evaluated during one period. 

(iii) The limit S of the iterative process. Since the problem 
is non-linear, we use an iterative method and solve a given 
sequence of linear problems in each time step. We stop the 
process when the relative difference of two consecutive solu¬ 
tions is less than <5. 


We measure the uncertainty by the metric: 

max(u, Mtheory) = maX ju-’ - M^j,eory| i (d) 

where denotes the temperature in the j-th time step. 
Figure [CT] shows the dependence of this metric on the tetra¬ 
hedra volume Afl and the time step At; the limit <5 is a con¬ 
stant here. We checked not only the maximum difference 
of u — Mtheory but also its mean dispersion, nonetheless, the 
results were very similar. We see that the uncertainty in¬ 
deed decreases with a refinement of the discretization. We 
use P/At = 1000 and Afl”^ = 100000 in our analyses. 

For high values of the thermal conductivity K, the 
changes of temperature are small and the iterative proce¬ 
dure converges quickly — the limit 5 = 10“"^ was reached 
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after only two or three iterations. For lower values of K, the 
procednre with the same S took about 10 iterations. As we 
foreshadowed in Section the procedure did not converge 
for the lowest values of considered range of K and we had 
to solve the problem using the relaxation method. 

Moreover, we have two possibilities of a posteriori un¬ 
certainty estimates. First, we can compare our results to 
Golubov and Krugly (2012 I (see Appendix E . Although our 


and their approaches are substantially different, both results 
are indeed comparable. Second, the mean torque induced 
by a boulder must vanish in the limit of zero conductivity, 
A —>■ 0, as we mentioned in Section |3.3| This is especially 
important as the low-A case cannot be tested with the ana¬ 
lytical solution ( All| |. Luckily, low K effectively corresponds 
to larger boulders (l/L) which do not contribute much to 
the total torque, so we consider this case as being of lesser 
importance. 


This paper has been typeset from a IAI]bX file prepared 

by the author. 
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